tau=.1;    % time step

ggap=70;

T=100;
time=tau:tau:T;

type=[];
type(end+1)=-1;      % soma
type(end+(1:2))=1;   % C_1
n=length(type);      % the number of compartments

m=max(type)-min(type)+1;  % the number of compartment types

% fs cell

mvr=-55+zeros(1,m);
mvt=[-40,-40+zeros(1,m-1)];
mp=[1,1+zeros(1,m-1)];
mC=[20,20+zeros(1,m-1)];
ma=[0.15,0.15+zeros(1,m-1)];
mb=[8,8+zeros(1,m-1)];
mc=[-55,-55+zeros(1,m-1)];
md=[200,200+zeros(1,m-1)];
mspike_peak=[25,25+zeros(1,m-1)];



vr=mvr(2+type);
vt=mvt(2+type);
p=mp(2+type);
C=mC(2+type);
a=ma(2+type);
b=mb(2+type);
c=mc(2+type);
d=md(2+type);

spike_peak=mspike_peak(2+type);

mgapm=ggap*[1 1 1 1 1 1];      % conductance (sensitivity) to j-1 (mother) compartment
mgapd=ggap*[1 1 1 1 1 1];     % conductance (sensitivity) to j+1 (daughter) compartment
gapm=mgapm(2+type);
gapd=mgapd(2+type);


nbrs=zeros(n,3);
nbrs(1,:)=[1 2 3];     % self-gap junctions (no effect), C_0
nbrs(2,:)=[1 2 2];     % C_s, self (no effect), C_0
nbrs(3,:)=[1 3 3];     % C_0, C_1, C_1


v=zeros(length(time),n);
v(1,:)=vr;
u=0*v;
g_ampa=0+0*v;
I=0+0*v;
st=65/tau;

decay = exp(-tau/6);

compartments = [1 2 3];
   g_ampa(floor(1/tau), 2 )=4.0;
   g_ampa(floor(30/tau), 3 )=4.0;
   g_ampa(floor(60/tau), 2 )=4.0;
   g_ampa(floor(60/tau), 3 )=3.9;
   
%   C(1:17) = [100 85 70 55  55 40 40 40 40 25 25 25 25 25 25 25 25 ] ;
%   p(1:17) = [3   2.5 2 1.5 1.5 1 1  1  1  1   1  1  1  1  1  1  1 ] ;
   
   
   for i=2:length(time)
       g_ampa(i,compartments)=g_ampa(i-1,compartments)*decay+g_ampa(i,compartments);
%       g_ampa(i,25)=300/2*t*exp(-2*t);
   end;

for i=1:length(time)-1    
 dv=p.*(v(i,:)-vr).*(v(i,:)-vt) -u(i,:) + I(i,:) + g_ampa(i,:).*(0-v(i,:)) + ...
     gapm.*(v(i,nbrs(:,1))-v(i,:)) + ...
     gapd.*(v(i,nbrs(:,2))+v(i,nbrs(:,3))-2*v(i,:));
 v(i+1,:)=v(i,:)+tau*dv./C;
 u(i+1,:)=u(i,:)+tau*(a.*(b.*(v(i,:)-vr)-u(i,:)));
 fired =find(v(i+1,:)>spike_peak);
 v(i,fired)=spike_peak(fired);
 v(i+1,fired)=c(fired);
 u(i+1,fired)=u(i+1,fired)+d(fired);
end;

plot(v(:,compartments))


legend(num2str(compartments(1)),num2str(compartments(2)),num2str(compartments(3)))
